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We solve the mixing-demixing transition in repulsive one-dimensional bose-bose mixtures. This is 
done numerically by means of the continuous matrix product states variational ansatz. We show that 
the effective low-energy bosonization theory is able to detect the transition whenever the Luttinger 
parameters are exactly computed. We further characterize the transition by calculating the ground- 
state energy density, the field-field fluctuations and the density correlations. 
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I. INTRODUCTION 


Understanding the physics of strongly correlated 
many-body systems is a formidable task, both in the lat¬ 
tice and in the continuum [1]. There is a fruitful synergy 
between condensed matter, high energy physics or quan¬ 
tum chemistry and the quantum information community. 
Ideas as tensor network states [2] or quantum simulations 
[3] are pursuing the goal of understanding phases and dy¬ 
namics beyond the paradigm of perturbative theories. 

One-dimensional many-body systems are a good ex¬ 
ample of this cooperation. Well established theoretical 
techniques as bosonization [4] are complemented with the 
density matrix renormalization group (DMRG) and ma¬ 
trix product states (MPS) [5, 6] in the lattice and more 
recently, continuous matrix product states (cMPS) in the 
continuum [7]. Ultracold gases are a paradigmatic ex¬ 
ample of experiments realizing one dimensional quantum 
fields [8, 9, 19]. Experiments and simulations in one di¬ 
mension are perfect testbeds since each of them can be 
used for benchmarking the other [10-13]. 

Of special relevance for this work is the cMPS formal¬ 
ism. Introduced by Verstraete and Cirac [7], these states 
constitute a variational class for the efficient simulation 
of quantum held theories that does not rely on a space 
discretization. The ansatz has proven to be efficient for 
computing the ground state, dispersion relation [14] and 
quantum evolution [15] of nonrelativistic theories. In ad¬ 
dition, introducing a suitable regularization prescription, 
it has also been applied to the study of certain relativis¬ 
tic phenomena [16]. Finally, the cMPS ansatz has been 
already tested for bose-bose [17] and fermi-fermi [18] mix¬ 
tures, both in nonrelativistic setups. 

Among the different scenarios covered in experiments, 
we are interested here in bose-bose mixtures [20, 21] . In 
this work, we aim to characterize the mixing/demixing 
phase transition occurring in repulsive bosonic mixtures, 
via cMPS. Here, the competition between the intra and 
interspecies couplings leads to the formation of two differ¬ 
ent phases. The miscible, were the two gases coexist, and 
the immiscible, where the two gases separate from each 


other [22]. This transition has been studied analytically 
within different approaches [23-25] but its numerical sim¬ 
ulation has been elusive. 

We will review the phase transition, clarify some issues 
on the instability occurring already in the Luttinger liq¬ 
uid description, and fully characterize it via the ground 
state energy, fluctuations and correlation functions. 


II. THE PHASE TRANSITION 


Two ID bosonic gases interacting via a quartic contact 
potential {h = 2m = 1) are described in second quanti¬ 
zation via the Hamiltonian 



where (^^(a:)) are the bosonic field operators 

which create (annihilate) bosonic particles of species a 
at the position x G [0, L]. They satisfy the commu¬ 
tation relations: ['0Q,(a;),'0^(a;')] = 6apd{x — x') and 

[ij)a{x),ipi3{x')\ = [-!/]]],(a:),'!^^(a:')] = 0. In this work, we 
want to characterize numerically the mixing-demixing 
transition occurring in mixtures whenever two differ¬ 
ent repulsive bosonic species are trapped together. We 
will consider the case where the participating species 
are nonconvertible, i.e. the individual particle densi¬ 
ties poa of each bosonic species are conserved sepa¬ 
rately [poa = {'4’tc{x)'^a{x)), { ) means averages over the 
ground state of (1)]. We will also restrict to the sym¬ 
metric case poi = Po 2 = Po, gii =922 = c > 0 and 
912 = 921 = fl /2 > 0 . 

The mixing/demixing transition has been broadly 
studied analytically, see e.g. Refs. 23-28. The phase 
separation, which lies on the competition between the re¬ 
pulsion strengths c and g, can be understood on several 
grounds. The simplest approach considers a mean-field 
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treatment. Here, the interaction term can be seen as a 
quadratic form of the densities. The latter is positive 
defined as long as g < 2c. When the positivity condi¬ 
tion is violated {g > 2c) an instability occurs. In more 
than one dimension, both species must separate in order 
to make the overlap integral zero, i.e., minimizing the 
repulsive interaction which yields the instability. In this 
phase {g > 2c) both species are immiscible as observed by 
resolving the spatial density profiles in the trap [22]. In 
ID, there is no possibility of spatial separation. The in¬ 
terspecies fluctuations, ipl{x)ipl{x)tpi{x)ip 2 ix) in Eq. (1) 
are zero in the ground state. 

One-dimensional systems are somehow special. Their 
confinement provides an enhancement of the collective 
behaviour, leading to a universality class in the low- 
energy or long wavelength sector. The latter is known 
as Luttinger Liquid [4, 29, 30]. This regime is de¬ 
scribed by introducing the bosonic operators 4>a and 9a 
in terms of which we rewrite the field operators '4’a {x) = 
(poc-5,0a(x)/^)i/2 Y.t=-oo This 

is nothing but the harmonic fluid approach treatment 
best known in the literature as bosonization [31]. Note 
that for high enough values of p, the exponential terms 
oscillate very fast and rapidly average to zero. Therefore, 
in order to obtain the low-energy effective Hamiltonian, 
we should only keep a few relevant terms. This leads to 


2nHes 


+VaKa{dJaf^ ( 2 ) 


stress that this result deviates from the mean field value 
g* = 2c. 

The phase separation has also been studied analyti¬ 
cally beyond perturbation theory by Kolezhuk [25]. He 
found that for one and two-dimensional gases, the tran¬ 
sition point, in the symmetric case, does not depend 
on the particle densities. Surprisingly enough, the non- 
perturbative result coincides with the mean-field descrip¬ 
tion. That is, the two species demix when g > g* = 2c. 
Following this result, one might be tempted to think that 
the bosonization framework is not able to predict cor¬ 
rectly the transition point. It could be argued that the 
Luttinger liquid paradigm breaks down at intermediates 
values of g, below the critical value g* = 2c. Here, we 
will show that this is not the case. We demonstrate that 
the bosonization predicts the transition correctly when 
the Luttinger parameters are computed exactly instead 
of using approximations. 


III. CMPS SOLUTION 


A translational invariant cMPS of N species {N = 2 
in this work) of bosonic particles is defined by the state 
vector [16]: 

dx Q (g) I -I- ^ (g) f’tixn 1^2) 

a=l ) 

(4) 


lx) = Ti'auxT’exp f J 


+ 



‘ 2 -gxdx^idxf 2 + gc cos (2{(j)i - ^ 2 ) 


This long wavelength description is fully characterized 
by the dimensionless parameters ATq,, the velocities Va 
and the coupling strengths gx and gc (Luttinger pa¬ 
rameters). For the symmetric case considered here, we 
have that Vi = V 2 = v and Ki = K 2 = K. This 
model can be easily decoupled by introducing the nor¬ 
mal modes = l/-\/2((^i ±(^ 2 ) and 9± = l/-\/2(0i ±02)- 
In terms of them, the low-energy Hamiltonian reads 

27riLeff = J dx J2u=± +v^K^{dx9^)'^'j + 

gcj dx cos(-\/8(^_). The normal modes’ velocities z;± are 
defined as follows 


4 = 1 ± 


Kgx 

V 


(3) 


As pointed out by Cazalilla and Ho in Ref. 24, the cou¬ 
pled system (2) is unstable when vf becomes negative. In 
other words, the action is not anymore definite positive, 
pretty much like in the mean-field argument sketched be¬ 
fore. This will happen whenever Kgx > v. Thus, to 
compute the transition point, we just need to find the 
Luttinger parameters from the original Hamiltonian (1). 
In Ref. 24, AT, gx and v were approximated via expres¬ 
sions valid in the weak interspecies coupling g regime. In 
the quasi-condensate regime ^ = cfp < 1 [31], the insta¬ 
bility is estimated to happen at g* = 2c(l — ^/7/27r). We 


where ipaix) are the bosonic field operators, Q and Ra 
are a set of complex, D x D matrices acting on an aux¬ 
iliary D-dimensional space and jH) is the free vacuum 
state vector (^a(x)|H) — 0). V denotes a path-ordering 
prescription and the partial trace, Traux, is taken over 
the auxiliary space. This way of writing field states is 
the continuous limit of a MPS [7]. Here, the dimension 
D of the auxiliary matrices corresponds to the so-called 
bond dimension, an upper bound to the entanglement 
entropy. Typically, the low-energy states of local Hamil¬ 
tonians should possess a low amount of entanglement, 
consequently D is a small number. Being the bond di¬ 
mension small, the state (4) represents an efficient trial 
for finding the ground state of one-dimensional field the¬ 
ories numerically. 

In a previous work [17] the authors showed how to 
construct a two-species cMPS starting from two decou¬ 
pled single species solutions. In brief, for coupled fields 
we considered coupled auxiliary spaces (one per bosonic 
field). The total auxiliary Hamiltonian was extended to 

(K = -tQ - ^ZRlRa) 

P 

K = Ki®l2+li®K2 + Y. ® ■ ( 5 ) 

p=i 

where Ka is the auxiliary Hamiltonian associated to 
bosonic species a. The parameter P accounts for the 
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FIG. 1. (Color online) Instability in the bosonization descrip¬ 
tion. The square velocities ni, defined by (7), are calculated 
using cMPS for different values of the parameter 7 = c/p: 0.52 
(filled circles), 1.5 (filled triangles), 2.38 (open squares) and 
3.0 (open diamonds). In the inset, we compare the numeri¬ 
cal result for 7 = 2.38 (open squares) with a weak-coupling 
estimation (dashed line) for the same value of 7 with c = 1.5 
and p = 0.63 (see the main text). All of the simulations have 
been performed with D = 5 and P = 1. 


number of pairs of coupling matrices entering in the 
cMPS state. Consequently, the matrices Ra, belonging 
to the auxiliary space of field a, were extended into the 
total product space: Ri — Ri I and i ?2 = I G) i? 2 - 
Denoting D the dimension of the matrices i?i and i ?2 
the bond dimension is then D = D^. The total number 
of variatonal parameters is D^(4 -b 2P). Details can be 
found in Ref. 17. 

In the thermodynamic limit (L —> 00 ), the fluctuations 
and correlation functions can be computed from 

Capix -y) = (-!/l(a:)V:^(j/)' 0 p(?/)Vi„(a:)) ( 6 ) 

= Tr[(Rp 0 R})e'^^--y\R^ ® R*)] , 

without loss of generality, we have assumed that x > y. 
Remind that, throughout this work ( ) means averages 
over the ground state of (1). The transfer operator T is 
defined as: T = Q(SjI + I(SjQ* Ra^^ Ra- Finally, 

note that the fluctuations are calculated by making x = y 
in ( 6 ). 


IV. RESULTS 


As already anticipated, our goal is to characterize the 
mixing/demixing transition numerically. We do it in two 
ways. First, we study the instability in the low-energy 
regime described by the effective Hamiltonian (2). The 
second strategy is to look directly at the ground state of 
( 1 ) and compute the fluctuations and correlation func¬ 
tions, Cf. Eq. ( 6 ). 


A. Bosonization instability 


In the harmonic fluid approach the normal modes for 
the fields decouple (see the discussion below Eq. (2)). 
Each of these modes propagate with different velocities, 
v±. Within the bosonization framework, these velocities 
can be related to the ground state energy density (cq). 
The explicit expressions for the velocities are [31, 32] 


= 2 p± 


d'^eo 

dpi 


(7) 


with p± = Pi ± p 2 - Analytical estimations for these ve¬ 
locities follow from Eq. (3). In the weak-coupling regime 
(g ^ c), it is safe to assume that v and K correspond to 
the solutions for a single bosonic field [31]. In turn, is 
approximated by g^ — gj'x (see reference 24). 

In the inset of figure 1, it can be seen that already 
at intermediate values of g (well below the critical value 
p*) the predicted velocities v± using weak-coupling an¬ 
alytical expressions deviate from the numerically com¬ 
puted ones [17, 32]. A consequence of this deviation is 
the failure on the estimation of the point where v2 be¬ 
comes negative, which in turns marks the critical value 
p». In the main figure 1 we have zoomed the v2 around 
the transition point for different values of 7 = c/p. 
As it has been already pointed out, within the weak- 
coupling treatment, the transition is estimated to hap¬ 
pen at p* = 2c(l — ^/7/27r). As 7 = c/p, the latter result 
makes the transition point dependent on both, the in¬ 
traspecies coupling c and the particle density p. On the 
other hand, once is exactly derived from the ground- 
state energy density by using relation (7), we see how the 
transition point becomes independent of 7 . In fact, we 
see that the mode propagating with velocity V- becomes 
ill-defined at g^fjc = 2 [33], in agreement with the mean- 
field and Kolezhuk results [25]. Therefore, once the Lut- 
tinger parameters are exactly computed, the bosoniza¬ 
tion predicts correctly the transition. 


B. Characterization beyond bosonization 

Having a full knowledge of the ground state of (1), we 
proceed now to characterize the phase transition beyond 
the bosonization formalism. In figure 2, we see the be¬ 
haviour of the ground state energy density as a function 
of the interspecies coupling. It is direct to realize that af¬ 
ter p*/c = 2 the energy remains constant. In this region, 
the ground state is such that the last term of ( 1 ), be., the 
one accounting for the interaction among different fields, 
has zero average. In other words, after the transition 
we have that: Ci, 2 ( 0 ) = (i/][(a;)Vji(a:)V>|(a;)V[ 2 (a:)) = 0 , 
which is explicitly represented in the inset of figure 2 
(open squares). This confirms our previous exposition 
for the phase transition: in one dimension, phase separa¬ 
tion implies zero interspecies fluctuations. 

Apart from the transition point estimation and the 
zero field-field overlapping nature for the demixed phase. 
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FIG. 2. (Color online) Ground state energy density of Hamil¬ 
tonian (1) as a function of the coupling ratio gjc calculated 
with cMPS. We keep fixed the intraspecies coupling c = 1.5 
while the particle density in each of the gases is equal to 
p = 0.5. In the inset we show the ground state fluctuations 
as a function of p/c: (711(0) (open triangles), (722(0) (inverted 
filled triangles) and (7i2(0) (open squares) defined in (6). Sim¬ 
ulations have been performed with D = 5 and P = 1. 

we can go further in characterizing the properties of the 
ground state before and after the transition. Let us start 
with the mixed phase. By looking at the inset of hg- 
ure 2 we see that the fluctuations ( 71 , 2 ( 0 ) do not re¬ 
main constant as soon as the interaction is switched on. 
The latter behaviour reflects a sublinear growth of eo in 
terms of g. This means that a simple mean field the¬ 
ory (■!/>J(a:)'(/)i(x)^|(x)'(/' 2 (a:^)) = P 1 P 2 is not sufficient for 
describing this phase. 

We will discuss now the demixed phase. As explained 
above, after the transition ( 7 i, 2 ( 0 ) = 0. It is straightfor¬ 
ward to see that a ground state of the form 

|7^dm) = ^ (lX 2 p) O \^) + O |x 2 p)) ( 8 ) 

fulfils this condition (suffix dm stands for demixing). Be¬ 
sides, iTfdin) niust satisfy the particle density conserva¬ 
tion for each bosonic species: {Xdni\i’t{^)i’a{x)\Xdni) = 
p, which in turn imposes that: (x 2 p|%(a^)V'a((c)|X 2 p) = 
2p. Indeed, this is confirmed in Fig. 2 via the ground 
state energy density. We see that after the transition, 
eo is the energy of a single bosonic gas (Lieb-Liniger 
model) with self-interaction c but double particle density 
2p (dashed line) [33]. Finally, by looking at the fluctua¬ 
tions Caa{0), we check that they coincide with those of 
a single bosonic gas with self-interaction c and particle 
density 2p, divided by a factor of two due to normaliza¬ 
tion in (8) . The fluctuations of a single gas are shown in 
the inset of figure 2 with a dashed line. 

We hnish our phase characterization by studying the 
correlation functions, Capix). The results are plotted 
in figure 3. By definition, the correlations at zero dis¬ 
tance match the fluctuations. On the other hand, in 
the limit x —>■ 00 , the correlations factorize yielding 



FIG. 3. (Golor online) Density correlation functions: Cii(x) 
(open triangles), C 22 (x) (inverted filled triangles) and Ci 2 (x) 
(open squares) as a function of the distance x for the same 
parameters of Fig. 2. The transition happens at p*/c = 2. 
We plot the correlations (a) before the transition g/c = 0.52 
and (b) after the transition g/c = 2.53. The shape of this 
curve brings to mind the popular story of the boa constrictor 
digesting an elephant [34]. Simulations have been performed 
with D = 5 and P = 1. 


Cap(x —> 00 ) = [marked with dashed lines in 3 a) 
and b)]. In the mixed phase the correlation length is of 
the order of a: = 5, pretty much the same than for a 
single bosonic species with self-interaction c and particle 
density p. 

More structure for Ca,p{x) appears in the demixed 
phase. The interspecies correlation function Ci 2 (x), ob¬ 
viously starting at zero, has a large correlation length 
~ 10"*^ (notice the logarithmic scale). To understand this 
large correlation length we recall that after the transition 
the fields are infinitely repelled. Our interpretation is re¬ 
inforced by looking at Caa{x). In the range 0 < x < 10 
the correlations build up to 2p^ wich means that they can 
be approximated by Caa{x) = l/ 2 (x 2 p|yOa|X 2 p)^ = 2 p^. 
Therefore, the coherence has been lost at the single field 
level. However, the fully uncorrelated state will involve 
the full state |7fdm) and pretty much like for the Ci 2 {x) 
correlations, the demixed phase is equivalent to an infi¬ 
nite repulsive phase, explaining again the large coherence 
lenght to reach the asymptotic limit Caa{x —>■ 00 ) = p^. 


V. CONCLUSIONS 

Summarizing, by means of cMPS we have computed 
numerically the ground state of two repulsive one¬ 
dimensional bosonic nonconvertible helds. This kind of 
systems exhibits the so-called mixing/demixing phase 
transition. We have validated previous analytical results 
for the transition point. Furthermore, we have demon¬ 
strated that this point can be resolved within the Lut- 
tinger liquid formalism whenever the effective parame¬ 
ters of the theory are calculated exactly. All this marks 
a step forward for the cMPS method, here, resolving a 
phase transition in a non-trivial quantum field theory. 
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